



DEPARTMENT OF MECHANICAL ENGINEERING & MECHANICS 
COLLEGE OF ENGINEERING & TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


7 ^ 25 * 


NUMERICAL METHODS FOR INCOMPRESSIBLE VISCOUS 
FLOWS WITH ENGINEERING APPLICATIONS 


By 

M. E. Rose, Co-Principal Investigator 
and 

R. L. Ash, Principal Investigator 
Final Report 

For the period ended April 16, 1986 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Research Grant NAS1-17993 
Task Authorization No. 32 

Mr. John B. Anders, Technical Monitor 
HSAD-Viscous Flow Branch 


CO 

o 

in 

CT> 

t— 

m o 

<N 

(0 i o 

1 

r-t M3 

CO 

u m 

CO 

a 

25 

a o 


ion 

CL 20D G3/34 

o 

» G '77 

Cm 

-O -H U 

CO 

u a 
o o 

Q 

CUQ 

O 

a o> 

SC 

h « 

fH 

H H 

W 

3rlO 

a 


t-i 

w a 
a 


O Cm M3 

u 

H CO 

H 

Cm CO O' 

03 

2 T— 

W 

in o 

a 

a H • 

a 

O H U 

a 

u a 


m u *c 
H M 
> h3 M3 

CO 

CM r- 


W CM 

*— 

■< O’- 

CO 

po a a. 

CO 

M O'rl 


m a ns m 

1 

in h am 

03 

WW (H 

u 

os W 

1 

a wfl ^ 

*£ 

a z o • 

in 

o H -H f> 

■=S 

UO U’H 

a 

2 2 CD R 


H w o. a 





August 1988 


DEPARTMENT OF MECHANICAL ENGINEERING & MECHANICS 
COLLEGE OF ENGINEERING & TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


NUMERICAL METHODS FOR INCOMPRESSIBLE VISCOUS 
FLOWS WITH ENGINEERING APPLICATIONS 

By 

M. E. Rose, Co-Principal Investigator 
and 

R. L. Ash, Principal Investigator 


Final Report 

For the period ended April 16, 1986 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Research Grant NAS1-17993 
Task Authorization No. 32 

Mr. John B. Anders, Technical Monitor 
HSAD-Viscous Flow Branch 


Submitted by the 

Old Doninion University Research Foundation 
P. 0. Box 6369 
Norfolk, Virginia 23508 


August 1988 
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FLOWS WITH ENGINEERING APPLICATIONS 


By 

M. E. Rose 1 and R. L. Ash 2 
ABSTRACT 

A numerical scheme has been developed to solve the incompressible, 
three-dimensional Navier-Stokes equations using velocity-vorticity vari- 
ables. This report summarizes the development of the numerical approxima- 
tion schemes for the divergence and curl of the velocity vector fields and 
the development of compact schemes for handling boundary and initial-boun- 
dary value problems. 
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2 Eminent Professor, Department of Mechanical Engineering & Mechanics, Old 
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INTRODUCTION 


This report is intended to document the research which was conducted in 
support of the development, a velocity-vorticity, Navier-Stokes solver 
(Ref.l). The nimerical techniques developed in this effort have been used 
subsequently in the study of three-dimensional vortex breakdown (Ref. 2). 

Approximation Schemes for div £ = 0, curl jj = £ 

3 

Let the planes x_ = const, describe a Cartesian grid in R , and denote 
by e^ a volume element with center at x^ . Suppose that the fundamental 
domain, D, is the union of such elements. As Fig. 1 shows, each element, e, 
has faces, y, edges, a, and vertices, v, and we identify these by a point 
associated with each. In addition, |ej, |y | , |a| denote the respective 

o 

volume, area, and length such that, if |o | = 0(h), then |y| = 0(h ), and |e| 
= 0(h 3 ). 

Consistent with this geometric construction, let u^e), u(y), and ^(a) 
denote the average values of jj on e, y, and a, respectively. The values 
ujv ) associated with the vertices are sometimes called box-variables and are 
often useful for quadrature evaluations of ]j(e), jj( y ) , and j^(o). For 
simplicity of notation, these will also designate the quadrature evaluation 
of these average values in terms of box-variables. 

Referring to Figure 1, consider the faces of an element e. . where 

1 , J , K 

1 2 
the index i is associated with the x or x-axis, j is associated with the x 

3 

or y-axis, and k is associated with the x or z-axis. The defining 
relations for the average and difference operators on the y.^^ j k ^ aces 
are 


^i,j,k> * = ^ (Y i + l/2,j,k ) + ^ (Y i-l/2,j,k^ /2 


( 1 ) 


f 



0 


Figure 1. Cartesian element e, with faces y, edges a, and vertices v. 


2 





( 2 ) 


‘l^i.O.k 1 ’^Hl/Z.J.k* - " (Y i-l/2J,k>l /2; 


the operators 2 anc * W 3’ A 3 are sim ^ ar ^ defined. The operators 
<S.u(y. . . ) are determined by 

1 — 1 , J , K 

b- = h 1 6 i , i = 1,2,3 (3) 

where h 1 = — — . The defining relations for the average and difference 
2 

operators relating sides and edges, and edges and vertices are similar. 

With these operator definitions in mind, we can construct approximation 
schemes for both div jj = 0 and curl _u = £ as follows: 

First, define the volume average of div _u over an element e as 


div„ u = / div u de. 

e - lel 


(4) 


Gauss's theorem evaluates div g u in terms of _u(y ) • jn on 3e, where in is the 
unit outward normal, i.e.. 


div u = u • n d |y | 

1 e | 9e 

= T~ § ii(y ) * 1 (5) 

I e J ye3e 

where y is the oriented area. By suitably arranging the order of summation 
and using Eqs. (2) and (3), Eq. (5) can be written as 


3 



Using box-variables to evaluate the average values on the right hand side of 
Eq. (6), one obtains 


div u 
e — 


+ « 2 u lU 3 u 2 + « 3 p 1 u 2 u 3 ). 


(7) 


In a similar fashion, define the surface average of the normal 
component of curl u over a surface y as 


n • curl u = / n • curl u dlyl* 

Y “ |y|y y “ 


( 8 ) 


Stokes' theorem evaluates n • 
the unit tangent vector, i.e.. 


curl . 
Y 


u in terms of u(a) 


a on 3 , where a is 
~ Y - 


n • curl^ u^ 


— — § u(a ) 

|y| ae3y 


a* 


(9) 


For a Cartesian element the opposite edges a + on a face y have equal 
lengths, and the associated unit tangent vectors have different signs. This 
allows one to rewrite the summation in Eq. (9) in terms of an operator 
acting on the edges. For example, the component on the x^-face is given by 


(n* curl^ £)| 1 = 6 2 m 3 (o 3 ) - 6 3 u 2 (a)« 


( 10 ) 


Using box variables to evalulate the average values on the right hand side 
of Eq. (10), one obtains 



(11) 


(n • cur 1 ^j) | x = ( 3 2 y 3 M 3 ” a 3 w 2 u 2^* 

With these defining relations, the general form for the (normal) vorticity 
components on the faces of an element, e, can be written as 

(n.* cur 1^ u)\ l = 6 j k u k ) “ 5 k^j u j^’ (12) 

where (i,j,k) indicates an even permutation of (1,2,3). 

Using the definitions (7) and (12), Fix and Rose (Ref. 3) have shown 
that the Cauchy-Rieman type equations described in (Ref. 1) can be solved by 
the least-squares solution of 

div u = 0 in D 
e — 

curl u = 8 on r 
r - 

£ • u. = n • u f 
when div 0 £ =‘ 0. 

Compact Schemes for Boundary and Initial-Boundary Value Problems 

The following discussion will outline the general development of the 
schemes that were used in developing the nunerical schemes for (Ref. 1). 
These schemes will be seen to provide domain decomposition extensions of 
conventional boundary integral and boundary element methods. As a result, 
both boundary value and initial-boundary value problems can be handled. 


( 12 ) 

(13a) 

(13b) 
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Consider a vector function, V, having s components. A typical feature 
of boundary value problems for systems of elliptic equations involving V on 
a domain, D, is that s-s' components of V on 3D are determined, by means of 
the. solution operator on D, by s' prescribed components of V on 3D, s' < s. 
We may call the s' prescribed components the primary variables and the s-s 1 
components the complementary variables . In certain simple cases the rela- 
tionship between these variables can be described by means of simple inte- 
gral equations. 

By constructing an approximate solution operator on D it may be possi- 
ble to determine the relationship between the primary and complementary 
variables at N points of interpolation on 3D. We may expect that these 
approximate values will converge, as N + «, to the solution values V under 
reasonable precautions about the construction. This is the basis of di- 
screte boundary integral methods. 

Returning to the continuous problem on D, suppose D is partitioned into 
volume elements, D = {e}. With arbitrary values of the primary variables 
chosen on the boundary of each element (but consistent with the values pre- 
scribed on 3e 3D) we can solve the boundary value problem in each element; 
the solution will be identical to the solution values of the boundary value 
problem in D in corresponding volume elements if they both have the same 
values on the boundary of each element and, thus, is continuous across in- 
terelement boundaries. 

This suggests the following discrete approximation method: In each 

element, e, choose the center point of each face, y, of e as an interpola- 
tion point and, using an appropriate solution operator in e, obtain the 
discrete boundary integral relationship between the primary and complemen- 
tary variables at the interpolation points on the boundary of 3e. We call 
this a compact equation on the element. Next impose continuity conditions 


6 



at the interpolation points in D and use values prescribed by the problem on 
D when the interpolation point lies on 3D. Then, solve the resulting 
algebraic problem. This, in effect, provides a domain decomposition 
extension of boundary element methods. (In this construction one may 
incorporate the continuity conditions quite simply by identifying the left 
and right limits at an interpolation point by their common value.) 

We call this construction a compact scheme . It requires for its de- 
velopment only an element-by-element description of the discrete integral 
equations which relate the primary and complementary variables at the inter- 
polation points on the boundary of each element e. As shown below, this 
idea may be applied to time dependent problems as well, (cf. Refs. 4-6). 

The weak-element method (Ref. 7) implements this construction by using 
as an approximation basis a manifold of solutions of the differential 
equation (or an approximation to it) in each element. Necessarily, then, 
the compact scheme which results is consistent with the differential 
equation in each element. This construction also leads immediately to a 
discrete energy estimate which approximates that which applies to the 
differential equation on D. Thus the convergence of the scheme is assured 
and leads to second order accurate results. 

We will now indicate how a simple Galerkin method can be used to obtain 
compact schemes for general volume elements. 

A Boundary Value Problem 

2 

As an example, we will discuss the Poisson equation v v = f. Consis- 
tent with earlier notations, let f(e) indicate the value at the center of e 
and v (y ) the value at the interpolation point on a face y. Define the 
bilinear boundary operator B g (v,w) by 
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so that Green's theorem can be written 


2 2 

B (v,w) = $ (vtf v - w w) d | e 1 « 


(15) 


If w is any solution of the homogeneous problem (v w = 0), then 


B p (v,w) = / wf d | e J = (w,f) 


(16) 


Second order accurate quadrature approximations to Eqs. 14 and 15 yield 


B^(v,w) = § [w( Y ) 

e L 


ye3 1 


3 n 


v ( y ) 3 -— - - 1 | y | 

3n -I 


(17) 


(w,f) = w(e)f (e) |e |* 


(18) 


Suppose e has i faces: let w, , i = 1,2 ,...,£ denote, say, the first i 

^ / 

2 

harmonic polynomials (i.e., V = 0). Compact equations on e are given 

by 


Bg (v,w.) = w is f)g, i = 1,2, ... ,£ 


(19) 


and provide an 0(h ; truncation error. 


These equations establish an algebraic relationship between the 1 



values v(y) and 9 - v — - on the faces of 3e. The coefficients in the equation 


an 


I » i / 

are determined by evaluating w. (y)|y| and [y| on the faces, i = 

1 an 

2 

1,2,. ..,4 . For the Poisson equation v v = f, application of these ideas is 
more straightforward if it is written as the system. 


v • p = f 


(20a) 


p = VV. 


(20b) 


seen, in the case of a Cartesian grid, use of the functions 


2 2 2 2 . 
w = (1, x, y, z, x - y , x - z ) 


( 21 ) 


in Eq. (19) leads to the compact equations 


6 p + 6 p +6 p =f 
x x y r y z z 


(22a) 


up = 6 V, up = 6 V, up = <5 V 
xx x y y y z z z 


(22b, c,d) 


u x v “ 7 h x 6 x p x = V ' f h ?yPy = y z v “ 7 h z s Z P Z 


1 ■- 2 - 

2 


1 ' 2 - 

2 


(22e,f) 


In order to show the convergence of this scheme, it is possible to 
construct a discrete energy estimate. Recall that the weak-element method 
(Ref. 7) constructs an approximate solution v on each element as 


9 



(23) 


e 

v = 


k 


c • v • 
n v i 


V, 


^ 2 

where visa particular solution of v v = f(e). Equation (19) determines 

6 

the coefficients by a Galerkin construction on each element. Since v is a 
2 e 

solution of V v = f(e) on e, it satisfies the energy equation 


* v e iL. d | y | - (v e , f) + / |vv e | 2 d|e|* 
3e 3 n e 


(24) 


Using the second order accurate quadrature formulas to evaluate the integral 
terms, one can approximate Eq. (24) as 

0 

I v e (r) ii-kl |y| = (» 2 , f) + ; (|vv e | 2 ) d|e| (25) 

ye3e 3n e 

Recalling the continuity conditions imposed by the compact construction, and 
summing over elements in D, we obtain 


0 

§ v e (Y) iXili |y| ■ I |{v e , f) ♦ / (|vv e | 2 ) d|e|). (26) 

ye3D 3n eeD e 


This is a discrete approximation to energy estimates for the solution in D, 
viz. , 


$ v !X d|y | = / (vf + |vv| 2 ) d | e | • (27) 

3D 3n. D 
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This discrete energy estimate, together with the obvious fact that the 
approximation v is consistent with the differential equation, implies by 
standard arguments that the scheme converges and, in fact, with second order 
accuracy. 

It is, of course, possible to extend these ideas to more general 

boundary value problems of the form Lv = f. For sufficiently small volume 

elements, L can be approximated by an operator with constant coefficients in 

each element, which we call L . In this case Green's theorem can be written 

e 

as (cf. Eq. (15)). 

B g (v,w) = / (wL g v - vL g w) d | e | (28) 

e 

* * 
where L g is the adjoint of l_ e . Solutions of the adjoint equation L w = 0, 

are easily generated in the form 

w = exp [a_ • x], (29) 

where a_ satisfies the characteristic polynomial equation L e (cO = 0. The 
compact equations which result, 

B (v,w) = (w, f) (30) 

e e 

may now involve exponential factors. It is possible to avoid the use of 

* 

exponentials by using polynomial solutions of L w = 0 which can be generated 
by 


11 



Wi = — : — (exp[a_ • x])| aa Q, i = 0,1,2,...; j = 1,2,3. (31) 

9 o 1 (j) 

The weak-element construction just described is based upon using a 
projection on a manifold of solutions of the differential equation in each 
element. This same idea can be applied to time-dependent problems as well. 

Initial-Boundary Value Problems 

As an example, we will base our discussion on the diffusion equation 
v t =v2v » e D, 0 < t < T (32a) 

with initial and boundary conditions 

v(x, 0) = g(x), t = 0, (32b) 


v()^, t) = v r (x r ), Xj, e 3D. 


(32c) 


From the discussion in (Ref. 1), it is sufficient to solve this problem in a 
time strip S m : |t - tj < t; i.e., with initial data vU, anc * 


boundary data v on 3D x S m , we seek to determine 


3v(x r , t) 


on 3D x S m as 


3n 


well as v(x_, t ffl + 1/2 ), x e D. 

Introducing a domain decomposition D = {e}, we consider the same type 
of problem on each cylinder set e x S m . The general solution can be written 
as 


v(x^, t; o_) = / A(a_) exp(a_ • _x - 0t) ch , 


(33) 
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where 8 satisfies the dispersion relation 


8 = !a| 2 , (34) 

and A(a_) is determined from the initial and boundary conditions in e x S m . 

If e has £ faces, one can seek an approximate solution which interpolates to 
the initial and boundary conditions on the element in the form 

£ 

v(x, t) = l A(a , ) exp(ct. • x - 0 -t)« (35) 

i=0 -i-i 

Once again, the relationship between the primary and complementary 
variables for the discrete problem can be determined by a Galerkin procedure 
using an appropriate form of Green's theorem. Let w indicate a solution of 
the adjoint equation 

L*(w) = Wj. + V 2 w = 0. (36) 

The application of Green's theorem to this problem leads to the relation 
(cf. Eq. (30)). 

— (w, v) e = B e (v, w) • (37) 

dt 

The approximation which results by interpolation is then (cf. Eq. 19)). 

— (w,v)jj = b£(v,w). (38) 

dt 


13 



m 


A time average of Eq. (38) on S produces the equation 


5 t (w, v)£ = Bg(v m , /) 


where v m , w m indicate time averages over the strip S m and 


V ( V = (v(t m +T) - v(t m -t))/4t 


M t v (t m ^ = + x) + T ^ /2 - 


m 


m 


( 39 ) 


(40a) 

(40b) 


Once again, the choice of s' + 1 solutions v^ of the adjoint Eq. (36) will 
then determine the s' + 1 complementary solution values at points of 
interpolation in e x S m in terms of the s' + 1 primary solution values. 

The result is a set of compact equations for the problem. Write Eq. (32a) 
as the first order system 




(41) 


In this case, e. is the interval |x - x.j < h. The discrete mixed initial 

value problem on e x S m can be stated as: given v(e, t m _^ 2 ) as 

m 3v m (x r ) 

data and v (x^), x^ e 3 e, as boundary data, determine (i) , x f e 3e, 

3 X 

and (n) v(e, t^). 
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Three elementary polynomial solutions of the adjoint Eq. (36) are 

2 


Wi - ( w o ,w i» w 2) - ( t+ — — ) , 

(42) 

where the origin is taken at the center of e x S m . 
then leads to 

The compact Eq. (39) 

V (e ’ t m ) = B e (vm > w o) /4x 

(43a) 

0 = Bg(v m ,wJ l ) /AX 

(43b) 

u t v(e,t m ) = B e( vm * w 2) /Ax * 

(43c) 

These simplify to 


V = V 

(44a) 

V = V 

(44b) 

h^ 

and p t v = u x v - — 6 x p, 

(44c) 


where we have suppressed^the reference to e x S m . The compact scheme 
results by requiring that v and p be continuous across endpoints of the 
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intervals interior to D, using prescribed values of v on 3D. (Note that the 
space operators in Eq. (44) apply to the face values of (v,p) on the 
cylinder e x S m , while the time operators apply to the values of v on the 
upper and lower bases.) 

It is possible to obtain an energy estimate for the system described by 
Equations (44). Multiply Eq. (44a) by y^v and use both Eqs. (44b, c); the 
resulting equation is 


1 

2 


v 2 + (3 x v) 2 



= « x ( vp) . 


(45) 


Summing over the elements in D,'the discrete energy estimate for the 
approximation is 

2 

I V 2 AX + n (6 x v) 2 + JL (6 x p) 2 ]ax = vp| 6D . (46) 


This expression corresponds to that for the continuous problem 



2 dt D 


v 2 d x + / (v x ) 2 dx = vp| aD . 


(47) 


The compact equations (44) are obviously consistent with the 
differential equation and, in view of the discrete energy estimate Eq. (46), 
the Lax equivalence theorem implies the convergence of the compact scheme. 
This same argument holds in three dimensions using an arbitrary partition of 
the domain D, D = { E} . 
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Finally, it is straightforward to extend these ideas to the advection- 
diffusion equation 

2 

v t + V • (av) = vV v, (48) 


where _a and v are coefficients, which in the context of the physical 
problems of Ref. 1 are associated with the velocity and viscosity of a 
fluid, respectively. Once again, writing this as a first order system 
yields the equations 


v 


t 


= V 


• P 


(49a) 


p = vW - av 


(49b) 


The corresponding adjoint equations are. 


-w 


t 


= Vq + v 



(50a) 


q = vVw 


(50b) 


The resulting form of Green's theorem is now 


— ( w , v ) = B (v,w) = / ( wj) - vqj • ndy. (51) 

dt e e 3e 


If, in the time-strip S™, the coefficient a in each element e is frozen as 


17 



a^ = _a m , then the elementary solution of Eq. (50) is 

w = exp [a_ • x_ - 8 1] (52a) 

where 

3 = v |a | 2 + _a m • x (52b) 

and, again, the origin is taken as the centerpoint of e x S m . The discrete 
form of Eq. (50) is obtained by using appropriate interpolated values in the 
quadrature approximation to B g and taking the time average over S m . 

The values a = 0, a = -v”*a m in Eq. (51) lead to steady-state 
solutions, while x - a m t is a simple time-dependent polynomial solution. 

Thus the appropriate approximation basis for this simple advection-diffusion 


equation is 

m 

V (wV = (*» x_ a 1 " 1 * exp[ - -iJi ] ) . (53) 

The compact equations which result are 

6 t v = B e W 0^ /Ax (54a) 

-»j x a m u t v = Bg (v rn ,wJ V )/Ax (54b) 

[ (9 m )’^ sinh 0 m ]6 t v = Bg( v m ,w^) /ax (54c) 
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where 


6 m = a m h/v . (54d) 

When Eqs. (54a) and (54c) are combined (using (43)), (44c) can be replaced 
by 


0 = Bg(v m ,(e m )" 1 sinh e m - w!J) . 

Expanding the right hand sides of Eqs. (54a, b and d), one obtains the 
compact scheme 

6 v =6 p 
t x 

->\ a 'V = * , x' > - u v 

% P - h p « x P = V8 X V - U x 3 m (c x u x v - h o 8 x v) 

where 

.1 m , m N -l 
p = coth e - (e ) 


(55) 


(56a) 

(56b) 

(56c) 


(5 6d ) 


and 


i . m, m*. - 1 
c x = 1 - 4 x a (“x a ) » 


(56e) 
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The coefficient p given by (55d) controls the weighting given to upwind 
terms in the compact scheme. It can be consistently approximated by 


P (0 ) = 9/3, |e| < 3 

= sgn 9 , |e I > 3 


(57) 


whose use allows us to employ an exponential-type scheme without having, in 
fact, to calculate exponential terms. 

The extension to three dimensions results by using the basis 

w i = (1, x - a^t, y - a 2 t, z - a 3 t, expl-a^x/v), 

exp(-a 2 y/v), exp(-a 3 z/v ) ) . (58) 
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